SUBROUTINE RADBCHEXCH(array, &
&  ix0, ix1, iy0, iy1, is0, is1, &
&  ixx0, ixx1, iyy0, iyy1, iss0, iss1, &
&  npnormbound)

! Note, that  >=2 MPI-processes in both x and y directions are assumed
! in this algorithm


use MPI_VARIABLES, only : &
& nxsh, nysh, nssh, &
& size_MPI, comm3d, &
& xlbound, xrbound, &
& ylbound, yrbound, &
& horisbound, &
& coords, dims, isperiodic

implicit none

integer(4), intent(in) :: ix0, ix1, iy0, iy1, is0, is1
integer(4), intent(in) :: ixx0, ixx1, iyy0, iyy1, iss0, iss1
integer(4), intent(in) :: npnormbound ! number of points to be exchanged in
                                      ! normal to boundary direction

real(8), intent(inout) :: array(ix0:ix1,iy0:iy1,is0:is1)

! Local variables
integer(4), parameter :: nodata = -999
integer(4), parameter :: neighborhood = 1
integer(4) :: i0, i1, j0, j1, k0, k1

real(8), allocatable :: work_arrsectget(:,:,:)

if (.not. horisbound) i0 = nodata

if (xlbound) then
  if (ylbound) then
    i0 = ixx1 + 1
    i1 = ixx1 + 1
    j0 = iyy0
    j1 = iyy0 + npnormbound
    k0 = iss0
    k1 = iss1
  elseif (yrbound) then
    i0 = ixx1 + 1
    i1 = ixx1 + 1
    j0 = iyy1 - npnormbound
    j1 = iyy1
    k0 = iss0
    k1 = iss1
  else
    i0 = ixx0
    i1 = ixx0 + npnormbound
    j0 = iyy0 - 1
    j1 = iyy0 - 1
    k0 = iss0
    k1 = iss1   
  endif
endif 

if (xrbound) then
  if (ylbound) then
    i0 = ixx0 - 1
    i1 = ixx0 - 1
    j0 = iyy0
    j1 = iyy0 + npnormbound
    k0 = iss0
    k1 = iss1
  elseif (yrbound) then
    i0 = ixx0 - 1
    i1 = ixx0 - 1
    j0 = iyy1 - npnormbound
    j1 = iyy1
    k0 = iss0
    k1 = iss1
  else
    i0 = ixx1 - npnormbound
    i1 = ixx1
    j0 = iyy0 - 1
    j1 = iyy0 - 1
    k0 = iss0
    k1 = iss1   
  endif
endif 


if ( .not. (xlbound .or. xrbound)) then
  if (ylbound) then
    i0 = ixx0 - 1
    i1 = ixx0 - 1
    j0 = iyy0
    j1 = iyy0 + npnormbound
    k0 = iss0
    k1 = iss1   
  elseif (yrbound) then
    i0 = ixx0 - 1
    i1 = ixx0 - 1
    j0 = iyy1 - npnormbound
    j1 = iyy1
    k0 = iss0
    k1 = iss1   
  else
!   Internal domain processes do not need any data
    i0 = ixx0
    i1 = ixx0
    j0 = iyy0
    j1 = iyy0
    k0 = iss0
    k1 = iss0
  endif
endif

!if (xlbound .and. ylbound) then
!  write(*,*) 'array before 1 ', array(i0:i1,j0:j1,k0:k1)
!endif

! Getting the first section
allocate (work_arrsectget(i0:i1,j0:j1,k0:k1))
call ARRSECTGET &
& (array, work_arrsectget, &
&  ix0, ix1, iy0, iy1, is0, is1, &
&  ixx0, ixx1, iyy0, iyy1, iss0, iss1, &
&  i0, i1, j0, j1, k0, k1, &
&  nxsh, nysh, nssh, &
&  comm3d, size_MPI, coords, dims, &
&  neighborhood, isperiodic)
deallocate (work_arrsectget)

!if (xlbound .and. ylbound) then
!  write(*,*) 'array after 1 ', array(i0:i1,j0:j1,k0:k1)
!endif        

if (xlbound) then
  if (ylbound) then
    i0 = ixx0
    i1 = ixx0 + npnormbound
    j0 = iyy1 + 1 
    j1 = iyy1 + 1
    k0 = iss0
    k1 = iss1
  elseif (yrbound) then
    i0 = ixx0
    i1 = ixx0 + npnormbound
    j0 = iyy0 - 1
    j1 = iyy0 - 1
    k0 = iss0
    k1 = iss1   
  else
    i0 = ixx0
    i1 = ixx0 + npnormbound
    j0 = iyy1 + 1
    j1 = iyy1 + 1
    k0 = iss0
    k1 = iss1
  endif
endif

if (xrbound) then
  if (ylbound) then
    i0 = ixx1 - npnormbound
    i1 = ixx1
    j0 = iyy1 + 1 
    j1 = iyy1 + 1
    k0 = iss0
    k1 = iss1
  elseif (yrbound) then
    i0 = ixx1 - npnormbound
    i1 = ixx1
    j0 = iyy0 - 1
    j1 = iyy0 - 1
    k0 = iss0
    k1 = iss1   
  else
    i0 = ixx1 - npnormbound
    i1 = ixx1
    j0 = iyy1 + 1
    j1 = iyy1 + 1
    k0 = iss0
    k1 = iss1
  endif
endif

if ( .not. (xlbound .or. xrbound)) then
  if (ylbound) then
    i0 = ixx1 + 1
    i1 = ixx1 + 1
    j0 = iyy0
    j1 = iyy0 + npnormbound
    k0 = iss0
    k1 = iss1   
  elseif (yrbound) then
    i0 = ixx1 + 1
    i1 = ixx1 + 1
    j0 = iyy1 - npnormbound
    j1 = iyy1
    k0 = iss0
    k1 = iss1
  else
!   Internal domain processes do not need any data
    i0 = ixx0
    i1 = ixx0
    j0 = iyy0
    j1 = iyy0
    k0 = iss0
    k1 = iss0
  endif
endif

!if (xlbound .and. ylbound) then
!  write(*,*) 'array before 2 ', array(i0:i1,j0:j1,k0:k1)
!endif

! Getting second section
allocate (work_arrsectget(i0:i1,j0:j1,k0:k1))
call ARRSECTGET &
& (array, work_arrsectget, &
&  ix0, ix1, iy0, iy1, is0, is1, &
&  ixx0, ixx1, iyy0, iyy1, iss0, iss1, &
&  i0, i1, j0, j1, k0, k1, &
&  nxsh, nysh, nssh, &
&  comm3d, size_MPI, coords, dims, &
&  neighborhood, isperiodic)
deallocate (work_arrsectget)
!if (xlbound .and. ylbound) then
!  write(*,*) 'array after 2 ', array(i0:i1,j0:j1,k0:k1)
!endif        

END SUBROUTINE RADBCHEXCH
